knitr::opts_chunk$set(fig.width=5, fig.height=5) 

rm(list=ls())


options(scipen=999)

pkgs <- c("dplyr", "ggplot2", "stargazer", "car", "tidyverse", "extrafont", "haven", "coefplot", "plotrix", "sensemakr")
# A function to load all above packages. Install if they have not been installed.
usePackage <- function(p){
  for (pkg in p){
    if (!is.element(pkg, installed.packages()[,1]))
      install.packages(pkg, dep = TRUE, repos = "https://cloud.r-project.org/")
    require(pkg, character.only = TRUE)
  }
}
usePackage(pkgs)


d1 <- read.csv("Study 3 Data.csv")


d1$Duplicate <- ifelse(d1$Q_RelevantIDDuplicate == "TRUE", 1, 0)
d1$Bot <- ifelse(d1$Q_RecaptchaScore == "NA", 1, 0 )
d1$KickedOut <- ifelse(d1$AttnChk_2 != "None" | d1$Sex == "", 1, 0)

d0 <- subset(d1, is.na(Duplicate) & d1$Bot == 0 & d1$KickedOut == 0)



d0$Republican <- 0
d0$Republican[d0$PartyID == "Republican" | d0$PartyID2 == "Closer to Republican Party"] <- 1
d0$Democrat <- 0
d0$Democrat[d0$PartyID == "Democrat" | d0$PartyID2 == "Closer to Democratic Party"] <- 1

d0$Independent <- as.numeric(ifelse(d0$Democrat == 1 | d0$Republican == 1, 0, 1))


#d <- d0
d <- subset(d0, d0$Independent == 0 )

d$FTDemVoters <- d$AffPolPre_1
d$FTRepVoters <- d$AffPolPre_2


d$FTOutPartyPeople <- ifelse(d$Republican == 1, d$FTDemVoters, d$FTRepVoters)
d$FTInPartyPeople <- ifelse(d$Republican == 1, d$FTRepVoters, d$FTDemVoters)

d$AffPolPeople <- d$FTInPartyPeople - d$FTOutPartyPeople



## Demographics

d$Age <- d$age


d$Black <- as.numeric(ifelse(d$Race == "Black or African-American", 1, 0))
d$White <- as.numeric(ifelse(d$Race == "White", 1, 0))
d$Hispanic <- as.numeric(ifelse(d$Race == "Hispanic or Latino", 1, 0))

summary(as.factor(d$Sex))
d$Female <- ifelse(d$Sex == "Female", 1, 0)

## -3 = Extremely Conservative; 3 = Extremely Liberal
#summary(as.factor(d$Ideology))


## Democrats
summary(as.factor(d$Democrat))

## Republicans
summary(as.factor(d$Republican))


## Education: 1 = LT HS; 2 = High school; 3 = Other post HS; 4 = Some college; 5 = AS; 6 = BA; 7 = MA; 8 = PhD
d$Education <- as.numeric(ifelse(d$education == -3105 | d$education == 10, NA_real_, d$education))
summary(as.factor(d$Education))

## Income (increasing), 1 = less than 15k; 24 = more than 250k

d$Income <- as.numeric(ifelse(d$hhi == -3105, NA_real_, d$hhi))

## Region: 1 = Northeast; 2 = Midwest; 3 = South; 4 = West

d$Region <- d$region


library(dplyr)
library(kableExtra)


# Calculate summary statistics
demographic_summary <- d %>%
  summarise(
    Mean_Age = round(mean(Age, na.rm = TRUE), 2),
    Proportion_Black = round(mean(Black, na.rm = TRUE), 2),
    Proportion_White = round(mean(White, na.rm = TRUE), 2),
    Proportion_Hispanic = round(mean(Hispanic, na.rm = TRUE), 2),
    Proportion_Female = round(mean(Female, na.rm = TRUE), 2),
    Mean_Education = round(mean(Education, na.rm = TRUE), 2),
    Mean_Income = round(mean(Income, na.rm = TRUE), 2),
    Proportion_Democrat = round(mean(Democrat, na.rm = TRUE), 2),
    Proportion_Republican = round(mean(Republican, na.rm = TRUE), 2)
  )

# Transpose the table
transposed_summary <- t(demographic_summary)
colnames(transposed_summary) <- c("Summary Statistics")

# Convert the transposed summary to a data frame and add row names
transposed_summary_df <- as.data.frame(transposed_summary)
transposed_summary_df$Variable <- c("Mean Age", "Proportion Black", "Proportion White", "Proportion Hispanic", "Proportion Female", "Mean Education", "Mean Income", "Proportion Democrat", "Proportion Republican")

# Generate LaTeX table
latex_table <- transposed_summary_df %>%
  select(Variable, `Summary Statistics`) %>%
  kable("latex", booktabs = TRUE, align = "c", col.names = c("Variable", "Summary Statistics")) %>%
  kable_styling(full_width = FALSE)

# Print the table
cat(latex_table)

d$RedistTreatment <- as.factor(ifelse(d$DistrictCondition == "DistrictInParty", 1, 0))
d$JudgesTreatment <- as.factor(ifelse(d$JudgesCondition == "JudgesInParty", 1, 0))
d$JudgesTreatment2 <- ifelse(d$JudgesTreatment == 1, "Yes", "No")
d$WardTreatment <- as.numeric(ifelse(d$WardCondition == "WardInParty", 1, 0))
d$BoardSizeTreatment <- as.numeric(ifelse(d$BoardSizeCondition == "BoardSizeInParty", 1, 0))



d$JudgesResponse0 <- ifelse(d$JudgesTreatment == 1, d$JudgesInParty, d$JudgesControl)
d$JudgesResponse <- case_when(
  d$JudgesResponse0 == "Strongly support" ~ 3,
  d$JudgesResponse0 == "Support" ~ 2,
  d$JudgesResponse0 == "Slightly support" ~ 1,
  d$JudgesResponse0 == "Neither support nor oppose" ~ 0,
  d$JudgesResponse0 == "Slightly oppose" ~ -1,
  d$JudgesResponse0 == "Oppose" ~ -2,
  d$JudgesResponse0 == "Strongly oppose" ~ -3,
  TRUE ~ NA_real_
)
d$JudgesResponse <- (d$JudgesResponse + 3) / 6
d$JudgesResponse <- as.numeric(d$JudgesResponse)

b1 <- lm(d$JudgesResponse ~ d$JudgesTreatment)


d$RedistResponse0 <- ifelse(d$RedistTreatment == 1, d$DistrictInParty, d$DistrictControl)
d$RedistResponse <- case_when(
  d$RedistResponse0 == "Strongly support" ~ 3,
  d$RedistResponse0 == "Support" ~ 2,
  d$RedistResponse0 == "Slightly support" ~ 1,
  d$RedistResponse0 == "Neither support nor oppose" ~ 0,
  d$RedistResponse0 == "Slightly oppose" ~ -1,
  d$RedistResponse0 == "Oppose" ~ -2,
  d$RedistResponse0 == "Strongly oppose" ~ -3,
  TRUE ~ NA_real_
)
d$RedistResponse <- (d$RedistResponse + 3) / 6
d$RedistResponse <- as.numeric(d$RedistResponse)

d$RedistTreatment2 <- ifelse(d$RedistTreatment == 1, "Yes", "No")

b2<-lm(d$RedistResponse ~ d$RedistTreatment  )


d$WardResponse0 <- ifelse(d$WardTreatment == 1, d$WardInParty, d$WardControl)
d$WardResponse <- case_when(
  d$WardResponse0 == "Strongly support" ~ 3,
  d$WardResponse0 == "Support" ~ 2,
  d$WardResponse0 == "Slighlty support" ~ 1,
  d$WardResponse0 == "Neither support nor oppose" ~ 0,
  d$WardResponse0 == "Slightly oppose" ~ -1,
  d$WardResponse0 == "Oppose" ~ -2,
  d$WardResponse0 == "Strongly oppose" ~ -3,
  TRUE ~ NA_real_
)
d$WardResponse <- (d$WardResponse + 3) / 6
d$WardResponse <- as.numeric(d$WardResponse)


b3<-lm(d$WardResponse ~ d$WardTreatment  )


d$BoardSizeResponse0 <- ifelse(d$BoardSizeTreatment == 1, d$BoardSizeInParty, d$BoardSizeControl)
d$BoardSizeResponse <- case_when(
  d$BoardSizeResponse0 == "Strongly support" ~ 3,
  d$BoardSizeResponse0 == "Support" ~ 2,
  d$BoardSizeResponse0 == "Slightly support" ~ 1,
  d$BoardSizeResponse0 == "Neither support nor oppose" ~ 0,
  d$BoardSizeResponse0 == "Slightly oppose" ~ -1,
  d$BoardSizeResponse0 == "Oppose" ~ -2,
  d$BoardSizeResponse0 == "Strongly oppose" ~ -3,
  TRUE ~ NA_real_
)
d$BoardSizeResponse <- (d$BoardSizeResponse + 3) / 6
d$BoardSizeResponse <- as.numeric(d$BoardSizeResponse)


b4 <-lm(d$BoardSizeResponse ~ d$BoardSizeTreatment  )


### Second order beliefs


d$SOJudgesDems <- case_when(
  d$JudgesDems == "All" ~ 3,
  d$JudgesDems == "Almost all" ~ 2,
  d$JudgesDems == "Most" ~ 1,
  d$JudgesDems == "About half" ~ 0,
  d$JudgesDems == "Less than half" ~ -1,
  d$JudgesDems == "Almost none" ~ -2,
  d$JudgesDems == "None" ~ -3,
  TRUE ~ NA_real_
)
d$SOJudgesDems <- (d$SOJudgesDems + 3) / 6
d$SOJudgesDems <- as.numeric(d$SOJudgesDems)



d$SOJudgesReps <- case_when(
  d$JudgesReps == "All" ~ 3,
  d$JudgesReps == "Almost all" ~ 2,
  d$JudgesReps == "Most" ~ 1,
  d$JudgesReps == "About half" ~ 0,
  d$JudgesReps == "Less than half" ~ -1,
  d$JudgesReps == "Almost none" ~ -2,
  d$JudgesReps == "None" ~ -3,
  TRUE ~ NA_real_
)
d$SOJudgesReps <- (d$SOJudgesReps + 3) / 6
d$SOJudgesReps <- as.numeric(d$SOJudgesReps)

d$SOJudgesInParty <- ifelse(d$Democrat == 1, d$SOJudgesDems, d$SOJudgesReps)
d$SOJudgesOutParty <- ifelse(d$Democrat == 1, d$SOJudgesReps, d$SOJudgesDems)


d$SORedistDems <- case_when(
  d$DistrictDems == "All" ~ 3,
  d$DistrictDems == "Almost all" ~ 2,
  d$DistrictDems == "Most" ~ 1,
  d$DistrictDems == "About half" ~ 0,
  d$DistrictDems == "Less than half" ~ -1,
  d$DistrictDems == "Almost none" ~ -2,
  d$DistrictDems == "None" ~ -3,
  TRUE ~ NA_real_
)
d$SORedistDems <- (d$SORedistDems + 3) / 6
d$SORedistDems <- as.numeric(d$SORedistDems)



d$SORedistReps <- case_when(
  d$DistrictReps == "All" ~ 3,
  d$DistrictReps == "Almost all" ~ 2,
  d$DistrictReps == "Most" ~ 1,
  d$DistrictReps == "About half" ~ 0,
  d$DistrictReps == "Less than half" ~ -1,
  d$DistrictReps == "Almost none" ~ -2,
  d$DistrictReps == "None" ~ -3,
  TRUE ~ NA_real_
)
d$SORedistReps <- (d$SORedistReps + 3) / 6
d$SORedistReps <- as.numeric(d$SORedistReps)



d$SORedistInParty <- ifelse(d$Democrat == 1, d$SORedistDems, d$SORedistReps)
d$SORedistOutParty <- ifelse(d$Democrat == 1, d$SORedistReps, d$SORedistDems)


d$SOWardDems <- case_when(
  d$WardDems == "All" ~ 3,
  d$WardDems == "Almost all" ~ 2,
  d$WardDems == "Most" ~ 1,
  d$WardDems == "About half" ~ 0,
  d$WardDems == "Less than half" ~ -1,
  d$WardDems == "Almost none" ~ -2,
  d$WardDems == "None" ~ -3,
  TRUE ~ NA_real_
)
d$SOWardDems <- (d$SOWardDems + 3) / 6
d$SOWardDems <- as.numeric(d$SOWardDems)



d$SOWardReps <- case_when(
  d$WardReps == "All" ~ 3,
  d$WardReps == "Almost all" ~ 2,
  d$WardReps == "Most" ~ 1,
  d$WardReps == "About half" ~ 0,
  d$WardReps == "Less than half" ~ -1,
  d$WardReps == "Almost none" ~ -2,
  d$WardReps == "None" ~ -3,
  TRUE ~ NA_real_
)
d$SOWardReps <- (d$SOWardReps + 3) / 6
d$SOWardReps <- as.numeric(d$SOWardReps)



d$SOWardInParty <- ifelse(d$Democrat == 1, d$SOWardDems, d$SOWardReps)
d$SOWardOutParty <- ifelse(d$Democrat == 1, d$SOWardReps, d$SOWardDems)


d$SOBoardSizeDems <- case_when(
  d$BoardSizeDems == "All" ~ 3,
  d$BoardSizeDems == "Almost all" ~ 2,
  d$BoardSizeDems == "Most" ~ 1,
  d$BoardSizeDems == "About half" ~ 0,
  d$BoardSizeDems == "Less than half" ~ -1,
  d$BoardSizeDems == "Almost none" ~ -2,
  d$BoardSizeDems == "None" ~ -3,
  TRUE ~ NA_real_
)
d$SOBoardSizeDems <- (d$SOBoardSizeDems + 3) / 6
d$SOBoardSizeDems <- as.numeric(d$SOBoardSizeDems)



d$SOBoardSizeReps <- case_when(
  d$BoardSizeReps == "All" ~ 3,
  d$BoardSizeReps == "Almost all" ~ 2,
  d$BoardSizeReps == "Most" ~ 1,
  d$BoardSizeReps == "About half" ~ 0,
  d$BoardSizeReps == "Less than half" ~ -1,
  d$BoardSizeReps == "Almost none" ~ -2,
  d$BoardSizeReps == "None" ~ -3,
  TRUE ~ NA_real_
)
d$SOBoardSizeReps <- (d$SOBoardSizeReps + 3) / 6
d$SOBoardSizeReps <- as.numeric(d$SOBoardSizeReps)



d$SOBoardSizeInParty <- ifelse(d$Democrat == 1, d$SOBoardSizeDems, d$SOWardReps)
d$SOBoardSizeOutParty <- ifelse(d$Democrat == 1, d$SOWardReps, d$SOBoardSizeDems)



## Response_Type variables 
d$JResponse <- d$JudgesResponse
d$RResponse <- d$RedistResponse
d$WResponse <- d$WardResponse
d$BResponse <- d$BoardSizeResponse

## Treatment Type variables
d$TreatJ <- d$JudgesTreatment
d$TreatR <- d$RedistTreatment
d$TreatW <- d$WardTreatment
d$TreatB <- d$BoardSizeTreatment

y <- dplyr::mutate(d, ID = row_number())


## DVs to common stem
y$DV_1 <- y$RResponse
y$DV_2 <- y$JResponse
y$DV_3 <- y$SOJudgesInParty
y$DV_4 <- y$SOJudgesOutParty
y$DV_5 <- y$SORedistInParty
y$DV_6 <- y$SORedistOutParty
y$DV_7 <- y$WResponse
y$DV_8 <- y$BResponse
y$DV_9 <- y$SOWardInParty
y$DV_10 <- y$SOWardOutParty
y$DV_11 <- y$SOBoardSizeInParty
y$DV_12 <- y$SOBoardSizeOutParty

## Reshape to long
longS <- reshape(y, direction = "long", 
                 varying = list(c("DV_1","DV_2","DV_3","DV_4","DV_5","DV_6","DV_7","DV_8","DV_9","DV_10","DV_11","DV_12")))
longS$DV <- longS$time
longS$response <- longS$DV_1

## Code treatment variable
longS$treat <- ifelse(longS$DV == 1, as.numeric(longS$TreatR), 
                      ifelse(longS$DV == 2, as.numeric(longS$TreatJ), 
                             ifelse(longS$DV == 7, as.numeric(longS$TreatW),
                                    ifelse(longS$DV == 8, as.numeric(longS$TreatB), NA))))
longS$treat <- longS$treat - 1
table(longS$treat)

## Code in vs out party
longS$inorout <- ifelse(longS$DV %in% c(3,5,9,11), 1, 
                        ifelse(longS$DV %in% c(4,6,10,12), 0, NA))

## Code judge or redistrict or ward or boardsize
longS$JorRorWorB <- ifelse(longS$DV %in% c(3,4), 1, 
                           ifelse(longS$DV %in% c(5,6), 0, 
                                  ifelse(longS$DV %in% c(9,10), 2,
                                         ifelse(longS$DV %in% c(11,12), 3, NA))))

m_1 <- lm(response ~ treat*AffPolPeople + I(DV - 1), data = subset(longS, DV %in% c(1, 2, 7, 8)))
summary(m_1)

library(sandwich)
ses_1 <- sqrt(diag(vcovCL(m_1, cluster = longS$id[longS$DV %in% c(1:2, 7:8)])))
t_1 <- abs(coef(m_1) / ses_1)
round(cbind(coef(m_1), 
            ses_1,
            t_1,
            (1 - pt(t_1, m_1$df.residual)) * 2
),
3)

results_table <- round(cbind(coef(m_1), 
                             ses_1,
                             t_1,
                             (1 - pt(t_1, m_1$df.residual)) * 2
),
3)

# Convert to data frame and add row names
results_table_df <- as.data.frame(results_table)
rownames(results_table_df) <- c("(Intercept)", "treat", "I(DV - 1)")

# Generate LaTeX table with kable
latex_table <- kable(results_table_df, format = "latex", booktabs = TRUE,
                     col.names = c("Coefficient", "Clustered SE", "t-value", "p-value"))

# Print the table
cat(latex_table)
## Model for SO beliefs

m_2 <- lm(response ~ inorout + JorRorWorB, data = subset(longS, DV %in% c(3:6, 9:12)))
summary(m_2)

# clustered SEs
ses_2 <- sqrt(diag(vcovCL(m_2, cluster = longS$id[longS$DV %in% c(3:6, 9:12)])))
t_2 <- abs(coef(m_2) / ses_2)
round(cbind(coef(m_2), 
            ses_2,
            t_2,
            (1 - pt(t_2, m_2$df.residual))*2
),
3)


results_table_2 <- round(cbind(coef(m_2), 
                               ses_2,
                               t_2,
                               (1 - pt(t_2, m_2$df.residual)) * 2
),
3)

# Generate LaTeX table with stargazer
stargazer(results_table_2,
          title = "Regression Results",
          column.labels = c("Coefficient", "Clustered SE", "t-value", "p-value"),
          covariate.labels = c("(Intercept)", "inorout", "JorRorWorB"),
          summary = FALSE,
          type = "latex")



control_mean_1 <- coef(m_1)[1]
control_se_1 <- sqrt(diag(vcovCL(m_1, cluster = longS$id[longS$DV %in% c(1:2, 7:8)]))[1])

treat_mean_1 <- coef(m_1)[1] + coef(m_1)[2]
treat_se_1 <- sqrt(control_se_1^2 + (vcovCL(m_1, cluster = longS$id[longS$DV %in% c(1:2, 7:8)]))[1, 2]^2)


# Calculate mean and standard errors for intercept and "inorout" variable in m_2 model
control_mean_2 <- coef(m_2)[1]
control_se_2 <- sqrt(diag(vcovCL(m_2, cluster = longS$id[longS$DV %in% c(3:6, 9:12)]))[1])

inorout_mean_2 <- coef(m_2)[1] + coef(m_2)[2]
inorout_se_2 <- sqrt(control_se_2^2 + (vcovCL(m_2, cluster = longS$id[longS$DV %in% c(3:6, 9:12)]))[1, 2]^2)


# Create table
table_data <- data.frame(
  variable = c("Control", "Support the proposal", "Beliefs about in-party support", "Beliefs about out-party support"),
  mean = c(control_mean_1, treat_mean_1, inorout_mean_2, control_mean_2),
  se = c(control_se_1, treat_se_1, inorout_se_2, control_se_2)
)

mean_data3 <- as.data.frame(table_data)
mean_data3$variable <- factor(mean_data3$variable, levels = c("Control","Support the proposal", "Beliefs about in-party support",  "Beliefs about out-party support"))

custom_palette <- c("Control" = "#D55E00", "Support the proposal" = "#0072B2", "Beliefs about in-party support" = "#009E73", "Beliefs about out-party support" = "#CC79A7")

soplotgp <- ggplot(mean_data3, aes(x = variable, y = mean, fill = variable)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9)) +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black", position = "identity") +
  geom_errorbar(aes(ymin = mean -1.96 * se, ymax = mean + 1.96 * se),
                width = 0.2, position = position_dodge(width = 0.9)) +
  labs(title = "General population support for policy proposals\n that would reduce the power of the party in the minority",
       x = "", y = "Mean Response") +
  scale_x_discrete(labels = c("Control:\n Non-partisan Support", "Treatment:\n Partisan Support", "Perceived\n In-party Support", "Perceived\n Out-party Support")) +
  theme(axis.line.y = element_blank(), axis.ticks.y=element_blank()) +
  theme(text = element_text(size=10, family="LM Roman 10")) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(hjust = 0.5, family = "LM Roman 10")) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black")) +
  theme(legend.position = "none") +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1),  limits =  c(0, 1)) +
  scale_fill_manual(values = custom_palette) # Add custom color palette

soplotgp

ggsave(filename="soplotgp.jpeg", plot=soplotgp, device="jpeg", height=5.5, width=5.5, units="in", dpi=500)


## Partisan heterogeneity 

#longS$AffPol <- longS$AffPolPeople

#longS$AffPolS <- longS$AffPol/100

m_1 <- lm(response ~ treat*Republican  + I(DV - 1), data = subset(longS, DV %in% c(1, 2, 7, 8)))
summary(m_1)

#m_1 <- lm(response ~ treat + Republican + AffPolS  + I(DV - 1), data = subset(longS, DV %in% c(1, 2, 7, 8)))
#summary(m_1)

library(sandwich)
ses_1 <- sqrt(diag(vcovCL(m_1, cluster = longS$id[longS$DV %in% c(1:2, 7:8)])))
t_1 <- abs(coef(m_1) / ses_1)
round(cbind(coef(m_1), 
            ses_1,
            t_1,
            (1 - pt(t_1, m_1$df.residual)) * 2
),
3)

results_table <- round(cbind(coef(m_1), 
                             ses_1,
                             t_1,
                             (1 - pt(t_1, m_1$df.residual)) * 2
),
3)

# Generate LaTeX table with stargazer
stargazer(results_table,
          title = "Regression Results",
          column.labels = c("Coefficient", "Clustered SE", "t-value", "p-value"),
          covariate.labels = c("(Intercept)", "treat", "Republican", "I(DV - 1)"),
          summary = FALSE,
          type = "latex")

library(ggplot2)
library(dplyr)
library(broom)

# Calculate coefficients and their clustered standard errors
tidy_model_1 <- tidy(m_1, conf.int = TRUE)

# Calculate the clustered standard errors
cluster_se <- sqrt(diag(vcovCL(m_1, cluster = longS$id[longS$DV %in% c(1, 2, 7, 8)])))

# Calculate the t values
t_values <- abs(tidy_model_1$estimate / cluster_se)

# Calculate the p values
p_values <- 2 * (1 - pt(t_values, df = m_1$df.residual))

# Add the new standard errors, t values, and p values to the data frame
tidy_model_1 <- tidy_model_1 %>%
  mutate(`Clustered SE` = cluster_se,
         `t value` = t_values,
         `p value` = p_values)

# Rename variables for more intuitive plotting and filter out "(Intercept)" and "I(DV - 1)" coefficients
tidy_model_1 <- tidy_model_1 %>%
  rename(Coefficient = estimate, term = term) %>%
  filter(!term %in% c("(Intercept)", "I(DV - 1)")) %>%
  mutate(term = factor(term, levels = c("treat", "Republican", "AffPolS"),
                       labels = c("Treatment\n effect", "Republican\n identification", 
                                  "Affective\n polarization"))) %>%
  mutate(conf.low.90 = Coefficient - qnorm(0.95)*`Clustered SE`,
         conf.high.90 = Coefficient + qnorm(0.95)*`Clustered SE`,
         conf.low.95 = Coefficient - qnorm(0.975)*`Clustered SE`,
         conf.high.95 = Coefficient + qnorm(0.975)*`Clustered SE`)
# Generate the coefficient plot with 90% and 95% confidence intervals
cplot1 <- ggplot(tidy_model_1, aes(y = term, x = Coefficient)) +
  geom_vline(xintercept = 0, color = "black") +
  geom_point(color = "black") +
  geom_segment(aes(x = conf.low.90, xend = conf.high.90, y = term, yend = term), size = 1, color = "black") +
  geom_segment(aes(x = conf.low.95, xend = conf.high.95, y = term, yend = term), size = 0.5, color = "black") +
  labs(title = "Percent increase in\n support for policies\n subverting minority power", 
       x = "", 
       y = "") +
  scale_x_continuous(labels = scales::percent_format(accuracy = 2), limits = c(-.1, .3)) +
  theme_bw() +
  theme(text = element_text(size = 10, family = "LM Roman 10"),
        plot.title = element_text(hjust = 0.5, family = "LM Roman 10"),
        axis.line.y = element_blank(), 
        axis.ticks.y = element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"),
        legend.position = "none")

ggsave(filename="coefplotgp.jpeg", plot=cplot1, device="jpeg", height=5, width=5, units="in", dpi=500)

longS$outorin <-  ifelse(longS$inorout == 0, 1, 0)

m_2 <- lm(response ~ outorin + Republican + AffPolS, data = subset(longS, DV %in% c(3:6, 9:12)))
summary(m_2)

# clustered SEs
ses_2 <- sqrt(diag(vcovCL(m_2, cluster = longS$id[longS$DV %in% c(3:6, 9:12)])))
t_2 <- abs(coef(m_2) / ses_2)
round(cbind(coef(m_2), 
            ses_2,
            t_2,
            (1 - pt(t_2, m_2$df.residual))*2
),
3)

results_table_2 <- round(cbind(coef(m_2), 
                               ses_2,
                               t_2,
                               (1 - pt(t_2, m_2$df.residual)) * 2
),
3)

# Generate LaTeX table with stargazer
stargazer(results_table_2,
          title = "Regression Results",
          column.labels = c("Coefficient", "Clustered SE", "t-value", "p-value"),
          covariate.labels = c("(Intercept)", "outorin", "Republican", "AffPolS"),
          summary = FALSE,
          type = "latex")

tidy_model_2 <- tidy(m_2, conf.int = TRUE, conf.method = function(x) vcovCL(x, cluster = longS$id[longS$DV %in% c(3:6, 9:12)]))


# Rename variables for more intuitive plotting and filter out "(Intercept)" coefficient
tidy_model_2 <- tidy_model_2 %>%
  dplyr::rename(Coefficient = estimate, `Clustered SE` = std.error, term = term) %>%
  dplyr::filter(!term %in% c("(Intercept)")) %>%
  dplyr::mutate(term = factor(term, levels = c("outorin", "Republican", "AffPolS"),
                              labels = c("Asked about\n the out-party", "Republican\n identification", 
                                         "Affective\n polarization"))) %>%
  dplyr::mutate(conf.low.90 = Coefficient - qnorm(0.95)*`Clustered SE`,
                conf.high.90 = Coefficient + qnorm(0.95)*`Clustered SE`,
                conf.low.95 = Coefficient - qnorm(0.975)*`Clustered SE`,
                conf.high.95 = Coefficient + qnorm(0.975)*`Clustered SE`)

# Generate the coefficient plot with 90% and 95% confidence intervals
cplot2 <- ggplot(tidy_model_2, aes(y = term, x = Coefficient)) +
  geom_vline(xintercept = 0, color = "black") +
  geom_point(color = "black") +
  geom_segment(aes(x = conf.low.90, xend = conf.high.90, y = term, yend = term), size = 1, color = "black") +
  geom_segment(aes(x = conf.low.95, xend = conf.high.95, y = term, yend = term), size = 0.5, color = "black") +
  labs(title = "Percent increase in\n beliefs about who\n would support these policies", 
       x = "", 
       y = "") +
  scale_x_continuous(labels = scales::percent_format(accuracy = 2), limits = c(-.065, .3)) +
  theme_bw() +
  theme(text = element_text(size = 10, family = "LM Roman 10"),
        plot.title = element_text(hjust = 0.5, family = "LM Roman 10"),
        axis.line.y = element_blank(), 
        axis.ticks.y = element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"),
        legend.position = "none")
  # Calculate coefficients and their clustered standard errors



cplot2

ggsave(filename="coefplot2gp.jpeg", plot=cplot2, device="jpeg", height=5, width=5, units="in", dpi=500)


library(gridExtra)

# Set the width and height of each plot
width = unit(5, "cm")
height = unit(7, "cm")

# Combine the two plots side by side using grid.arrange()
cplotm <- grid.arrange(cplot1 + theme(legend.position = "none"), 
                       cplot2 + theme(legend.position = "none"), 
                       nrow = 1, 
                       widths = c(width, width), 
                       heights = c(height))

ggsave(filename="coefplotmgp.jpeg", plot=cplotm, device="jpeg", height=5.5, width=7, units="in", dpi=500)



### DEMOGRAPHICS MODEL


m_1 <- lm(response ~ treat + Republican + AffPolS + White + Black + Hispanic + Education + Income + Female + Age  + I(DV - 1), data = subset(longS, DV %in% c(1, 2, 7, 8)))
summary(m_1)

library(sandwich)
ses_1 <- sqrt(diag(vcovCL(m_1, cluster = longS$id[longS$DV %in% c(1:2, 7:8)])))
t_1 <- abs(coef(m_1) / ses_1)
round(cbind(coef(m_1), 
            ses_1,
            t_1,
            (1 - pt(t_1, m_1$df.residual)) * 2
),
3)

results_table <- round(cbind(coef(m_1), 
                             ses_1,
                             t_1,
                             (1 - pt(t_1, m_1$df.residual)) * 2
),
3)

# Generate LaTeX table with stargazer
stargazer(results_table,
          title = "Regression Results",
          column.labels = c("Coefficient", "Clustered SE", "t-value", "p-value"),
          covariate.labels = c("(Intercept)", "treat", "Republican", "AffPolS", "I(DV - 1)"),
          summary = FALSE,
          type = "latex")



### IDEOLOGY AND THREAT


longS$Ideologue <- ifelse(longS$Ideology == "Extremely conservative" | longS$Ideology == "Extremely Liberal", 1, 0)


m_1 <- lm(response ~ treat*Ideologue + I(DV - 1), data = subset(longS, DV %in% c(1, 2, 7, 8)))
summary(m_1)




